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Abstract 

A  transient  three-dimensional,  single-phase  and  non-isothermal  numerical  model  of  polymer  electrolyte  membrane  (PEM)  fuel  cell  with  high 
operating  temperature  has  been  developed  and  implemented  in  computational  fluid  dynamic  (CFD)  code.  The  model  accounts  for  transient 
convective  and  diffusive  transport,  and  allows  prediction  of  species  concentration.  Electrochemical  charge  double-layer  effect  is  considered.  Heat 
generation  according  to  electrochemical  reaction  and  ohmic  loss  are  involved.  Water  transportation  across  membrane  is  ignored  due  to  low  water 
electro-osmosis  drag  force  of  polymer  polybenzimidazole  (PBI)  membrane.  The  prediction  shows  transient  in  current  density  which  overshoots 
(undershoots)  the  stabilized  state  value  when  cell  voltage  is  abruptly  decreased  (increased).  The  result  shows  that  the  peak  of  overshoot  (undershoot) 
is  related  with  cathode  air  stoichiometric  mass  flow  rate  instead  of  anode  hydrogen  stoichiometric  mass  flow  rate.  Current  is  moved  smoothly  and 
there  are  no  overshoot  or  undershoot  with  the  influence  of  charge  double-layer  effect.  The  maximum  temperature  is  located  in  cathode  catalyst 
layer  and  both  fuel  cell  average  temperature  and  temperature  deviation  are  increased  with  increasing  of  current  load. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Polymer  electrolyte  membrane  (PEM)  fuel  cell  is  an  elec¬ 
trochemical  device  in  which  energy  of  chemical  reaction  is 
converted  into  electricity  directly.  Useful  features  such  as  high 
power  density,  simple,  safe  construction  and  fast  start-up  make 
those  particularly  suitable  for  home  appliance,  vehicles  and 
transportation  tools  [1].  Research  and  development  of  PEM  fuel 
cell  system  have  been  increased  dramatically  in  the  past  few 
years,  and  most  of  them  are  concerned  with  steady  state  con¬ 
dition.  However,  design  and  control  of  PEM  fuel  cell  will  also 
require  a  thorough  understanding  of  its  transient  characteristics 
which  are  important  for  power  conditioning  in  residential  appli¬ 
cations  and  for  automotive  systems.  Although  some  dynamic 
characteristics  can  be  determined  from  experiments,  a  model 
capable  of  predicting  transient  response  is  necessary  so  that 
system  behavior  can  be  analyzed  by  means  of  computer  simula¬ 
tion  in  different  conditions  of  load  current,  pressure  of  reactant 
gases,  temperature  and  stack  voltage,  etc.  Capability  of  pre¬ 
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dieting  transient  response  is  also  proved  useful  for  optimizing 
control  schemes  [2], 

Prior  to  this  work,  a  number  of  PEM  fuel  cell  transient 
mathematical  models  were  available  in  literature.  Almost  all 
these  models  were  based  on  system-level  which  were  simple 
and  included  all  important  characteristics  of  PEM  fuel  cell 
[3-9].  Cell  voltage  was  modeled  as  a  function  of  various  con¬ 
tributing  variables,  such  as  current,  cell  temperature  and  partial 
pressure  of  oxygen  at  cathode  inlet.  Effects  of  ohmic  resis¬ 
tance,  over  voltage  and  electrochemical  charge  double-layer 
capacitance  were  also  considered.  These  studies  established 
good  basis  for  understanding  transient  response  in  PEM  fuel 
cell.  However,  these  models  were  not  enough  to  understand 
intricate  actual  physical  processes  inside  PEM  fuel  cell  com¬ 
pletely  because  transient  behaviors  of  various  properties  were 
not  included  simultaneously,  particularly  the  governing  funda¬ 
mental  processes,  such  as  gas  transport  through  gas  diffusion 
layer  and  effect  of  electrical  conductivity.  Recently,  computa¬ 
tional  fluid  dynamics  (CFD)  and  improved  transport  models  had 
made  the  development  of  more  realistic  computational  mod¬ 
els  available  [10-16].  Detail  description  of  physical-chemical 
processes  and  coupled  transport  equations  with  electrochemi¬ 
cal  kinetics  were  included.  Most  of  these  models  were  focused 
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Nomenclature 

a  activity  of  gas  species 

/tjnem  cell  active  area  (m2) 

Cdf  differential  capacitance  of  porous  catalyst  layer 
(Fm-2) 

Cdi  effective  double-layer  capacitance  per  unit  geo¬ 
metric  volume  (F  m-3) 

Dkj  binary-diffusion  coefficient 

E  sum  of  internal  and  kinetic  energy  (J) 

£a  active  energy  (J) 

F  Faraday  constant  (96487  C  mol- 1 ) 

h  enthalpy  of  gas  species 

iref  reference  exchange  current  density  (A  m-3) 

!Qef  reference  exchange  current  density  at  reference 
temperature  (Am-3) 

/  local  current  density  (A  m-2) 

/aVg  average  current  density  (A  m-2) 

j  current  density  (A  m-3) 

J  diffusive  mass  flux  vector 

K  permeability  of  porous  electrode  (m-2) 

M  molecular  weight  of  gas  species  (kg  mol- 1 ) 
p  pressure  (Pa) 

R  universal  gas  constant  (8.3 14  J  mol- 1  K-1 ) 
e  entropy  (JK-1) 

Su  source  term  vector  of  momentum  equation 

Sm,  Sk,  >S'h,  S’soi ,  Smem  source  terms  of  continuity,  species, 
energy  and  charge  equation 
t  time  (s) 

T  temperature  (K) 

To  reference  temperature  (393  K) 

u  gas  flow  velocity  vector  (ms-1) 

Vo  thermodynamic  equilibrium  potential  (V) 

Vik  atomic  diffusion  volume 

Y  mass  fraction  of  gas  species 

Greek  symbols 

a  charge  transfer  coefficient 

P  empirically  determined  concentration  parameters 

s  porosity  of  the  electrode  materials 

<p  potential  (V) 

rj  local  over  potential  (V) 

k  electrical  (ionic)  conductivity  (S  m-1) 

ko  pre-exponential  factor  (S  m- 1 ) 

X  thermal  conductivity  ( W  m- 1  K- 1 ) 

Aeff  effective  thermal  conductivity  (W  m- 1  K- 1 ) 
p  dynamic  viscosity  (kg  s- 1  m-2) 

p  density  of  gas  mixture  (kg  m-3) 

r  viscous  stress  tensor 

reff  effective  stress  tensor 

§ss  specific  surface  area  of  porous  catalyst  layer 
(m3  m-2) 

pi  constant  number  determined  according  to  the 
experiment 


Subscripts  and  superscripts 
a  anode 

c  cathode 

f  fluid 

H2  hydrogen 

H2O  water 

K  species 

mem  membrane  electrolyte  phase 
O2  oxygen 

s  solid 

sol  electrode  solid  phase 


on  steady  state  behavior  in  order  to  facilitate  the  understand¬ 
ing  of  physical  process  and  mechanism  in  PEM  fuel  cell.  More 
recently,  Shimpalee  et  al.  [2]  and  Wang  and  Wang  [17]  extended 
the  single-domain  model  of  Um  et  al.  [18]  by  considering  all 
important  transient  processes  occurred  in  PEM  fuel  cell.  Tran¬ 
sient  response  of  cell  was  predicted  in  terms  of  local  distribution 
of  gas  concentration  and  current  density.  However,  effect  of 
electrochemical  charge  double-layer  was  ignored  according  to 
theoretical  estimation  of  various  time  constants  where  time  con¬ 
stant  of  charge  double-layer  was  sufficiently  short  to  be  ignored 
[17]. 

It  has  been  known  that  in  PEM  fuel  cell,  charge  double-layer 
is  occurred  in  a  thin  layer  on  or  near  electrochemical  reaction 
interface  and  acted  as  an  electrical  capacitor  during  transient 
procedure.  Such  capacitance  is  usually  in  the  order  of  a  few 
Farads,  which  is  high  in  terms  of  capacitance  value  [19].  This 
charge  double-layer  can  result  in  that  cell  current  is  moved  gently 
and  smoothly  to  new  value  in  response  to  step  change  of  cell 
operating  condition. 

In  this  paper,  a  complete  transient  numerical  model  capa¬ 
ble  of  characterizing  transient  response  of  high  temperature 
(T>  393  K)  PEM  fuel  cell  with  polymer  polybenzimidazole 
(PBI)  membrane  was  proposed  with  incorporating  capacitor 
effect  of  charge  double-layer.  Transient  characteristics  of  PEM 
fuel  cell  were  modeled  by  extending  the  three-dimensional 
model  in  previous  paper  [20]  with  accumulation  (unsteady) 
term.  Potential  field  equations  were  modified  with  considering 
about  charge  double-layer  capacitance  in  electrochemical  reac¬ 
tion  interface  based  on  macro-homogeneous  theory  of  porous 
electrode  [21,22],  Variation  of  charge  double-layer  capacitance 
with  potential  was  ignored.  All  the  important  transient  processes 
occurred  in  PEM  fuel  cell  with  high  operating  temperature,  such 
as  gas  transport,  species  dispersion  and  charge  double-layer  dis¬ 
charge  were  considered  comprehensively.  In  this  study  the  PEM 
fuel  cell  was  assumed  to  be  operated  at  433.15  K.  Product  water 
was  assumed  to  be  vaporous  and  treated  as  ideal  gas.  Water 
transportation  across  membrane  was  ignored  due  to  low  water 
electro-osmosis  drag  force  in  PBI  membrane.  The  model  was 
developed  and  implemented  in  the  commercial  CFD  code  FLU¬ 
ENT  6.1  with  custom  developed  user-defined-functions  (UDF) 
[23], 
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2.  Mathematical  model 


Fig.  1  schematically  shows  that  PEM  fuel  cell  can  be  divided 
into  following  sub-regions:  current  collector,  gas  flow  channel, 
gas  diffusion  layer  (GDL)  and  catalyst  layer  (CL)  in  both  anode 
and  cathode  sides,  and  membrane  in  the  middle.  Reactant  feed  is 
conveyed  by  gas  flow  channel  and  distributed  to  anode  and  cath¬ 
ode  side.  Then  reactant  is  passed  through  respective  porous  GDL 
and  reached  CL  where  electrochemical  reaction  is  occurred.  The 
membrane  acts  as  gas  separator,  electrolyte  and  proton  conduc¬ 
tor.  Electrons  are  collected  by  anode  current  collector  which 
is  connected  to  cathode  current  collector  through  external  cir¬ 
cuit.  In  this  study,  anode  feed  is  pure  dry  hydrogen  and  dry 
air  is  paralleled  in  cathode  gas  channel.  Since  gas  streams  in 
gas  flow  channel  are  at  low  velocity  (or  low  Reynolds  number), 


laminar  flow  with  ideal  gas  behavior  is  assumed.  Membrane  is 
considered  impermeable  to  gases.  Product  water  is  assumed  to  be 
vaporous  and  water  transportation  across  membrane  is  ignored. 

2.7.  Governing  equations 

Governing  equations  based  on  the  conservation  laws  of 
mass,  momentum,  species,  energy  and  charge  by  extending  the 
three-dimensional  model  of  [20]  with  including  time  dependent 
analysis  can  be  written,  in  vector  form,  as 

dp 

Continuity  e—  +  V  •  (pu)  =  Sm  (1) 

Momentum  -  +  -L y  .  (puu)  m  — V p  +  V  •  z  +  Sa 


Species  ediP*k)  +  V  -{puYk)  =  V  ■  Jk  +  Sk 


Energy  —  +  V  •  (u(pE  +  p)) 

=  V  •  ^AeffVJ  -  ^2 hkJk  +  (Teff  •  u)j  +  Sh 


Charge  -  Qi—  +  V  ■  (ksoiV0soi)  +  Ssoi  =  0 


■  0.6  and  0.4. 


(5) 
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Cdl  —  +  V  •  (/Cmem V 0me m)  +  5mem  =  0 
at 

p  is  density  of  gas  mixture,  which  is  defined  as 
1 

P  =  , 


(7) 


EkYk/Pk 

where  pk  is  density  of  species  k  and  it  can  be  calculated  according 
to  ideal  gas  law  relation 
pMk 


Pk  = 


RT 


(8) 


Jk  is  diffusive  mass  flux  vector,  which  can  be  written  as 
Jk  =  -Y.PDkjVYj  (9) 


Dkj  is  binary-diffusion  coefficient  which  can  be  calculated 
according  to  the  empirical  expression  [24] 


TlJ5(l/Mk  +  1/M;)1/2 

p((E,%)1'3  +  (E1t'e)1/s)2 


(10) 


where  Vik  is  atomic  diffusion  volume  and  the  value  of  V[k  is 
given  by  Cussler  [24],  Aeff  is  effective  thermal  conductivity  in 
a  porous  material  consisting  of  electrode  solid  matrix  and  gas, 
which  is  given  by 


A.eff  =  eAf  +  (1  —  s)Xs 


(ID 


where  as  is  thermal  conductivity  of  electrode  solid  matrix.  Af 
is  thermal  conductivity  of  gas,  which  can  be  expressed  as  a 
polynomial  function  of  temperature 

Xf  =  g0  +  SiT  +  g2T2  +  g3T3  (12) 

where  gi,  i  =  0, . . 3  can  be  determined  according  to  experi¬ 
ment.  Similar  to  the  thermal  conductivity,  the  viscosity  and  heat 
capacity  of  each  gas  species  can  also  be  described  by  polyno¬ 
mial  expression  of  temperature  with  empirical  coefficient.  More 
details  can  be  referred  in  [20]. 

Solid  electrode  phase  potential  dw  Eq.  (5),  which  is  solved  in 
current  collector,  GDL  and  CL  regions  of  both  anode  and  cath¬ 
ode  sides,  accounts  for  electron  transportation  through  electrode 
solid  conductive  materials.  Membrane  electrolyte  phase  poten¬ 
tial  0mem  Eq.  (6),  which  represents  proton  transportation  through 
CL  and  membrane,  is  solved  in  membrane  and  CL  regions. 
Charge  double-layer  effect  is  incorporated  with  time  dependent 
term.  The  double-layer  capacitance  per  unit  geometric  volume 
of  electrochemical  reaction  porous  electrode,  Cdi,  is  given  by 
[21] 


in  the  regions  except  CL.  p  is  local  over-potential  which  can  be 
expressed  as 

Pa  =  <psol  ~  0mem  (14) 

V c  =  d>sol  -  0mem  -  Vo  (15) 

Vo  is  the  thermodynamic  equilibrium  potential  which  can  be 
given  by  [19] 

V0  =  1.17  -  2.756  x  10“4(T  -  373.15) 

+  4.308  x  1(T5  In  ( m^a°2)  '  \  (16) 

y  oh2o  J 


This  equilibrium  potential  is  calculated  from  thermodynamic 
data  of  reaction  enthalpy  and  entropy  changes  while  product 
water  is  assumed  to  be  vaporous.  Definitions  of  an2,  a02  and 
«h2o  are 


«h2o 


P\ i2o 
Ph2o 


(17) 


where  pn2 ,  po2  and  ph2o  are  partial  pressure  of  gas  species. 

Detail  of  various  source  terms  in  Eqs.  (l)-(6)  can  be  seen 
in  Table  1.  It  shows  that  either  generation  or  consumption  of 
gas  species  and  creation  of  electric  current  are  occurred  only 
in  CL  region  where  electrochemical  reaction  is  taken  place. 
Current  densities  at  anode  and  cathode  can  be  described  by 
Butler- Volmer  equation  as  [26] 


exp 


(^) 


(18) 


.  =  ^»[«p(-^)-«p(^)]  «9, 

where  aa  and  ac  are  anodic  and  cathodic  charge  transfer  coeffi¬ 
cients,  see  in  Table  2.  i™f  and  are  reference  exchange  current 
densities  which  depend  on  local  temperature. 


'  =  L,o  exP  — 


Ea,  a  (  1 


'  =  L.o  exP  — 


£a,c 


(20) 

(21) 


where  E^  ,d  and  £a,c  are  active  energies  [22],  irael0  and  i™ q  are 
reference  exchange  current  densities  at  reference  temperature 
7o.  /3h2,  Po2  and  Pu2o  are  empirically  determined  concentra¬ 
tion  parameters  for  /1h2  =0.25,  fio2  =0.5  and  Ph2o  =  0.25. 
Local  current  density  I  in  membrane,  according  to  Ohm’s  law, 
is  defined  as 


Cdl  =  IssCdf 


(13) 


(22) 


where  £ss  is  specific  surface  area  of  CL  which  depends  on 
its  porous  structure,  and  Cdf  is  differential  capacitance  which 
depends  on  property  of  catalyst,  such  as  catalyst  loading  [25]. 
As  we  have  known  that  charge  double-layer  effect  is  occurred  at 
the  electrochemical  reaction  interface  between  electrolyte  and 
electrode.  In  this  study,  Cdi  is  assumed  to  be  homogeneous  in 
CL  of  both  anode  and  cathode  sides,  and  the  value  of  Cdi  is  zero 


Temperature  dependence  of  membrane  ionic  conductivity 
can  be  accurately  described  by  Arrhenius  equations  [27] 

=  (23) 

where  £a,k  is  activation  energy  and  kq  is  pre-exponential  factor 
at  reference  Tq.  Average  current  density,  which  is  the  average 
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Table  1 


Source  term; 

i  for  governing  equations  in  various  regions  of ; 

i  polymer  electrolyte  fuel  cell 

Source  terms 

Flow  channels 

GDL 

CL 

Membrane 

Continuity 

Sm,a  =  Sm,c=0 

Sm,a  =  Sm,c  =  0 

V  a  =  -^7a 

Sm,c  =  -  Tjr-  7c  +  -Jfh-  jc 

- 

Momentum 

Su,a=Su,c=0 

Su,a  =  SUjC  =  —gj^U 

5u>a  =  Su,c  = 

- 

Species 

Sfc,a  =  Sfcc=0 

Sk.i  =  Sk.c  =  o 

c  M«2  ■ 

So2'C  =  -^tjc 

Sh2o.c  =  -%£7c 

Energy 

Sh,a  =  Sh,c=0 

5'h,a  =  Sh,c  =  ® 

Sh,a  =  0 

Sh,c='-§T\As\  +  \jcVc\ 

5sol,a  =  Ssol,c  =  0 

Ssol,a=Ssol,c  =  ° 

Ssol,a  =  -7a 

.S’sol=0 

C”2' 

Smem,a=Smem,c=0 

Smem,a  =  Smem,c0 

5’sol.C  =jc 
£» m,a=7a 

W=-i 

Smem  =  0 

value  of  local  current  density  in  normal  direction  over  entire 
membrane,  can  be  obtained  by 

l avg  =  —  J  I  •  dA,  (24) 

where  Amem  is  cell  active  area. 

The  heat  released  from  CL  is  caused  by  enthalpy  change  and 
irreversibility  related  to  charge  transfer  [15].  In  this  paper,  ohmic 
heating  is  ignored  in  current  collector,  GDL  and  CL  but  consid¬ 
ered  in  membrane  due  to  its  relatively  low  ionic  conductivity. 
Empirically,  change  of  entropy  A.v,  as  a  function  of  temperature 
T,  can  be  expressed  as  [19] 

As  =  33.64  +  4.52564  x  \0~2T  -  2.98397  x  l(Tsr2 

+  3.40625  x  10“9T3  -  2.60417  x  KT12!4  (25) 

which  is  valid  for  temperature  from  373  to  1137  K. 


Table  2 

Electrochemical  properties 


Parameter 

Symbol 

Value 

Unit 

Porosity  of  GDL 

£ 

0.8 

- 

Porosity  of  CL 

£ 

0.6 

- 

GDL/CL  hydraulic  permeability 

5.0e-12 

m-2 

Membrane  ionic  conductivity 

*0 

12.99 

Sm-1 

GDL/CL  electric  conductivity 

K 

103.3 

Sm"1 

Current  collector  electric 

K 

535 

Sm"1 

conductivity 

Anodic  charge  transfer  coefficient 

a  a 

1.0 

- 

Cathodic  charge  transfer 

«c 

1.0 

- 

coefficient 

Anode  reference  exchange 

if0 

1.0e8 

Am"3 

current  density 

Cathode  reference  exchange 

*c0 

1.5e2 

Am-3 

current  density 

Thermal  conductivity  of  CL 

1.7 

W  m_1  K_1 

Thermal  conductivity  of 

0.95 

W  m_1  K_1 

membrane 

Thermal  conductivity  of  GDL 

1.7 

W  m_1  K_1 

Thermal  conductivity  of  current 

25 

W  m_1  K_1 

collector 

2.2.  Boundary  and  initial  conditions 

Eqs.  (l)-(6)  form  a  complete  set  of  governing  equations.  The 
corresponding  boundary  and  initial  conditions  were  described 
as  follows. 

(1)  On  gas  channel  inlet  boundaries,  stoichiometric  mass  flow 
rate  and  mass  fractions  of  species  were  prescribed  according 
to  the  current  with  cell  voltage  of  0.4  V.  On  outlet  bound¬ 
aries,  pressure  boundary  conditions  were  applied  with  cell 
operating  pressure  value  of  1.1  atm. 

(2)  In  Fig.  1,  on  both  anode  and  cathode  external  boundaries 
normal  to  y-axis,  which  contact  with  external  electrical  cir¬ 
cuit,  electrical  current  was  generated  only  through  these 
boundaries.  Step  change  of  cell  operating  voltage  between 
0.4  and  0.6  V  was  applied. 


Table  3 

Geometrical  and  operating  parameters 


Parameter 

Value 

Unit 

Cell  width 

3.4 

mm 

Channel  length 

235 

mm 

Channel  height 

0.7 

mm 

Anode  channel  width 

0.7 

mm 

Cathode  channel  width 

1.0 

mm 

Anode  GDL  thickness 

0.34 

mm 

Anode  CL  thickness 

0.04 

mm 

Membrane  thickness 

0.065 

mm 

Cathode  CL  thickness 

0.11 

mm 

Cathode  GDL  thickness 

0.34 

mm 

Thickness  of  current  collector 

2 

mm 

(both  anode  and  cathode) 

Operating  temperature 

433 

K 

Anode  stoichiometric  mass  flow  rate 

1.5 

- 

Cathode  stoichiometric  mass  flow  rate 

2.0 

- 

Anode  outlet  pressure 

1.1 

atm 

Cathode  outlet  pressure 

1.1 

atm 

Anode  inlet  mass  fraction  FL 

1.0 

- 

Cathode  inlet  mass  fraction  02:N2 

0.22:0.78 

- 

current  density  (A  cm’2) 
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(3)  Zero  flux  of  </>soi  was  defined  on  interface  between  CL  and 
membrane.  Zero  flux  of  was  defined  on  interface 
between  GDL  and  CL. 

(4)  In  isothermal  case,  cell  temperature  was  assumed  to  be 
constant  with  value  of  433.15  K.  In  non-isothermal  case, 
inlet  gas  was  assumed  to  be  pre-heated  with  temperature 
433. 15  K  and  outside  environment  temperature  was  also  set 
to  be  433.15  K. 

(5)  As  shown  in  Fig.  1,  the  fuel  cell  with  single  straight  gas 
channel  was  employed.  Generally,  fuel  cell  flow  plate  was 
formed  with  a  plurality  of  gas  channels.  Therefore,  symmet¬ 
rical  boundary  condition  was  applied  on  the  side  surfaces 


between  0.6  and  0.4  with  different  anode  and  cathode  stoichimometric  mass 


normal  to  x-axis.  No-slip  boundary  condition  was  applied 
on  remaining  walls. 

(6)  The  initial  conditions  were  applied  with  the  steady-state 
fields  for  cell  operating  voltage  of  0.6  V. 

2.3.  Numerical  procedure 

The  numerical  procedure  was  based  on  the  commercial  flow 
solver  FLUENT  6.1.  Volumetric  sources  or  sinks  in  Table  1  were 
implemented  in  the  solver  with  custom  developed  user-defined- 
function  (UDF).  Eqs.  (1)— (6)  were  solved  in  fully  implicit 
scheme  with  second-order  precision  for  both  spatial  and  tem¬ 
poral  discretization.  Therefore,  the  numerical  procedure  was 


between  0.6  and  0.4  with  effect  of  charge  double-layer. 


gas  flow  dir 


Fig.  8.  Local  transient  membrane  potential  contour  at  different  times  while  cell  voltage  changes  from  0.6  to  0.4  with  double-layer  considered  and  the  value  of  Cdi  is 
8.0  x  106Fm“3  (unit:  V). 


stable  throughout  the  time  step  At.  Note,  that  in  this  paper,  sev¬ 
eral  values  of  At  had  been  tested  in  order  to  make  sure  that 
the  results  were  independent  of  time  step  size.  Furthermore, 
stringent  numerical  tests  were  performed  during  steady  state 
simulations  to  ensure  that  the  solutions  were  independent  of 
grid  size.  A  mesh  with  about  180,000  grid  points  and  a  time  step 
size  At  of  0.005  s  were  found  to  provide  sufficient  spatial  and 
temporal  resolution. 

3.  Results  and  discussion 

As  shown  in  Fig.  1,  a  single-channel  PEM  fuel  cell  with 
polymer  polybenzimidazole  (PBI)  membrane  is  chosen  for  a 
parametric  study.  Values  of  relative  parameters  can  be  found  in 
Table  3. 

3.1.  Overshoot  and  undershoot 

Generally,  overshoot  or  undershoot  of  current  density  is  found 
during  step  change  in  cell  operating  condition,  such  as  cell  volt¬ 
age.  It  is  caused  by  the  fact  that  change  of  local  gas  concentration 
is  always  delayed  behind  change  of  local  potential. 

Fig.  2  shows  response  in  average  current  density  to  one  time 
step  change  in  cell  voltage  without  considering  about  the  influ¬ 
ence  of  charge  double-layer.  When  cell  voltage  is  changed  in 
one  time  from  0.6  to  0.4  V  at  f  =  0.0  s,  overshoot  is  occurred. 
The  average  current  density  is  rapidly  increased  from  0.467 
to  1.016  A  cm-2,  then  decreased  with  time  and  stabilized  with 
value  of  0.997  A  cm-2.  When  cell  voltage  is  changed  in  one  time 


from  0.4  back  to  0.6  V  at  t=  1.5  s,  undershoot  is  occurred  oppo¬ 
sitely.  The  average  current  density  is  rapidly  decreased  from 
0.997  to  0.448  A  cm-2,  then  increased  with  time  and  finally  sta¬ 
bilized  with  value  of  0.467  A  cm-2.  Fig.  3(a)  shows  local  current 
density  distribution  on  intermediate  cross-section  of  membrane 
at  four  different  time  steps  while  cell  voltage  is  changed  from 
0.6  to  0.4  V.  At  initial  sate  corresponding  to  cell  voltage  of  0.6  V 
at  time  f  =  0.0s,  the  local  current  density  near  gas  inlet  side  is 
high  and  decreased  along  gas  flow  direction.  Since  the  stoichio¬ 
metric  mass  flow  rate  of  anode  and  cathode  sides  are  prescribed 
as  1.25  and  2.0  when  cell  is  operated  at  voltage  of  0.4  V,  it  cor¬ 
responds  to  stoichiometric  mass  flow  rate  of  2.7  and  4.3  at  cell 
voltage  of  0.6  V.  We  have  large  excess  amount  of  fuel  and  air  at 
time  t  =  0.0  s.  Therefore,  the  current  density  distribution  across 
membrane  is  distinctively  different  from  those  (Fig.  3(a))  in  [20]. 
Variation  is  small  over  entire  cell  active  area.  The  maximum  and 
minimum  local  current  densities  are  0.490  and  0.413  A  cm-2, 
respectively.  When  cell  voltage  is  dropped  to  0.4  V  in  0.02  s,  the 
local  current  density  is  significantly  increased  over  entire  cell 
active  area  and  the  contour  pattern  is  distinctively  different  from 
that  at  t  =  0.0  s.  The  highest  local  current  density  is  1 .092  A  cm-2 
which  is  located  under  current  collector  land  area  near  gas  inlet 
side  because  of  the  influence  of  gas  concentration  and  ohmic 
losses  in  CL  and  GDL.  The  lowest  value  is  0.849  A  cm-2  and 
located  under  gas  channel  near  gas  outlet  side.  As  time  increases, 
local  current  density  is  decreased  while  the  distribution  contour 
pattern  is  similar  to  that  at  t  =  0.02  s.  At  time  t=  1 .0  s,  when  it  is 
stabilized,  the  highest  local  current  density  is  1.089  A  cm--  and 
the  lowest  value  is  0.818  A  cm-2. 
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between  0.6  and  0.4  with  non-isothermal  model. 

Fig.  3(b)  shows  transient  local  oxygen  mole  concentration 
contours  on  intermediate  cross-section  of  cathode  catalyst  layer. 
It  can  be  seen  that  local  oxygen  mole  concentration  distribution 
is  consistent  with  the  distribution  of  local  current  density.  Mole 
concentration  of  oxygen  is  decreased  monotonically  along  gas 
flow  direction.  The  concentration  value  under  current  collector 
land  area  is  smaller  than  that  under  gas  channel  due  to  gas  dif¬ 
fusion  and  consumption.  When  cell  voltage  is  dropped  to  0.4  V 
in  0.02  s,  it  can  be  seen  that  oxygen  is  depleted  rapidly  towards 
gas  outlet  side  and  the  concentration  distribution  becomes  much 
more  non-uniform.  As  time  increases,  the  oxygen  concentration 
is  decreased  gradually  especially  in  region  near  gas  outlet  side. 
Oxygen  mole  concentration  variation  over  entire  cathode  CL 
is  increased  from  1.50  mol  m-3  at  f  =  0.0s  to  2.95  mol  m-3  at 
t=  1.0  s.  The  oxygen  mole  concentration  near  gas  inlet  side  is 
stabilized  with  12.4%  reduction,  whereas  52.5%  reduction  near 
gas  outlet  side. 

When  cell  voltage  is  changed  from  0.4  back  to  0.6  V  at 
t  =  1.5  s,  on  the  contrary  to  Fig.  3(a),  local  current  density  is 
significantly  decreased  over  entire  cell  active  area  suddenly,  and 
then  increased.  Variation  of  local  current  density  over  entire 
cell  active  area  is  decreased.  This  is  caused  by  the  facts  that 
oxygen  distribution  variation  is  decreased  since  stoichiometric 
mass  flow  rate  is  increased  from  1.25  and  2.0  (at  0.4  V)  to  2.7 
and  4.3  (at  0.6  V). 

Fig.  4  shows  membrane  electrolyte  phase  potential  con¬ 
tours  at  midway  cross-section  area  of  membrane.  Due  to  the 
non-uniform  local  current  production  in  adjacent  CL,  there  are 
gradients  in  both  x  and  y  directions.  Since  charge  double-layer 
effect  is  not  considered  at  this  moment,  it  can  be  seen  that  elec¬ 
trolyte  phase  potential  distribution  can  be  stabilized  much  more 
promptly  than  oxygen  concentration. 

Fig.  5  presents  the  influence  of  gas  stoichimometric  mass 
flow  rate  on  cell  transient  response.  When  stoichiometric  mass 
flow  rate  of  air  is  increased  form  2.00  to  3.33  at  cell  voltage 
of  0.4  V,  average  current  density  is  increased  and  peaks  of  cur¬ 
rent  overshoot  and  undershoot  are  reduced  at  the  same  time.  In 
contrast,  when  stoichiometric  mass  flow  rate  of  air  is  decreased 
form  2.00  to  1 .25  at  cell  voltage  of  0.4  V,  average  current  density 
is  decreased  while  peaks  of  current  overshoot  and  undershoot 


are  enlarged.  Moreover,  form  the  result,  it  can  be  seen  that  sto¬ 
ichiometric  mass  flow  rate  of  hydrogen  does  not  influence  the 
cell  dynamic  performance.  This  is  caused  by  the  fact  that  anode 
side  is  fed  with  pure  dry  hydrogen.  It  is  distinctively  different 
from  [2],  where  the  water  activity  in  anode  side  is  considered. 

3.2.  Double-layer  effect 

Charge  double-layer  is  a  complex  electrode  phenomenon. 
In  electrochemical  systems,  a  layer  of  charge  on  or  near  elec¬ 
trode/electrolyte  interface  is  a  store  of  electrical  charge  and 
energy,  and  behaves  like  an  electrical  capacitor  [19],  In  Eq. 
(13),  the  double-layer  capacitance  per  unit  geometric  volume 
is  proportional  to  specific  surface  area  and  differential  capaci¬ 
tance.  The  specific  surface  area  f  ss  depends  on  porous  structure 
of  porous  CL  electrode  and  the  differential  capacitance  Qf 
depends  on  catalyst  content.  In  this  study,  value  of  §ss  is  set 
to  be  4125  m2  cm-3  [28]  and  value  of  Cdf  is  varied  from 
2x  10-2p,Fcm-2  and  2pFcm-2  [29],  Therefore,  according 
to  Eq.  (13),  value  of  Cdi  is  set  to  vary  from  8.25  x  105  Fm“3 
to  8.25  x  107  F  m-3.  In  this  paper,  property  of  both  anode  and 
cathode  CL  electrodes  are  assumed  to  be  identical.  This  is  not 
necessary,  and  one  might  identify  cases  where  non-identical 
electrodes  are  desirable. 

Fig.  6  shows  the  influence  of  charge  double-layer  on  cell 
transient  response.  It  can  be  seen  that  current  is  moved  smoothly 
with  increasing  value  of  Cdi.  From  prior  research  [17,21],  it 
is  well  known  that  time  constant  for  double-layer  charging  is 
proportional  to  value  of  Cdi.  Therefore,  it  takes  much  more  time 
for  the  cell  with  large  value  of  Cdi  to  stabilize  after  cell  voltage  is 
changed.  From  Fig.  6,  it  can  be  seen  that  value  of  current  density 
at  stabilized  state  is  independent  from  Cdi.  It  means  that  charge 
double-layer  just  influences  cell  transient  response  and  does  not 
influence  cell  steady  performance. 

Fig.  7(a)  shows  local  current  density  distribution  on  inter¬ 
mediate  cross-section  of  membrane  at  four  different  time  steps 
while  cell  voltage  is  changed  from  0.6  to  0.4  V.  Effect  of  charge 
double-layer  is  considered  with  Cdi  value  of  8.0  x  106  F  m-3. 
Since  charge  double-layer  just  influences  cell  transient  response, 
the  local  transient  current  density  contour  will  be  same  to  those  in 


voltage  changes  between  0.6  and  0.4. 
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Fig.  1 1 .  Temperature  distribution  contour  on  midway  section  at  different  times  while  cell  voltage  changes  from  0.6  to  0.4  without  considering  about  charge 
double-layer  (unit:  K). 


Fig.  3(a)  at  time  t  =  0.0  s  with  cell  voltage  0.6  V.  At  time  t  =  0.02  s 
and  cell  voltage  of  0.4  V,  comparing  with  Fig.  3,  due  to  influ¬ 
ence  of  charge  double-layer,  current  density  is  increased  more 
slowly  and  the  variation  is  smaller  than  those  without  consid¬ 
ering  about  charge  double-layer.  The  maximum  and  minimum 
local  current  densities  are  0.712  A  cm-2  and  0.522  A  cm-2.  As 
time  increases,  local  current  density  is  increased  smoothly  and 
the  value  near  gas  inlet  side  is  increased  faster  than  that  near 
gas  outlet  side.  In  Fig.  7(b),  as  local  current  density  is  increased 
smoothly  due  to  charge  double-layer  effect,  oxygen  mole  con¬ 
centration  is  decreased  also  smoothly  and  the  value  near  gas  inlet 
side  is  decreased  more  slowly  than  that  near  gas  outlet  side. 

Fig.  8  presents  membrane  electrolyte  phase  potential  con¬ 
tours  at  midway  cross-section  area  of  membrane  corresponding 
to  current  density  of  Fig.  7(a).  At  time  t  =  0.02  s,  distribution  con¬ 
tour  pattern  is  distinctively  different  from  that  in  Fig.  4,  which 
is  caused  by  charging  of  charge  double-layer.  As  time  increases, 
absolute  value  of  membrane  potential  is  increased  and  the  distri¬ 
bution  contour  pattern  is  become  similar  to  those  of  steady  state 
gradually. 

3.3.  Temperature  response 

In  previous  sections,  the  PEM  fuel  cell  is  assumed  to  be  oper¬ 
ated  at  constant  temperature  with  value  of 433. 15  K.  In  real  PEM 
fuel  cell,  heat  is  released  from  cathode  catalyst  layer  due  to  elec¬ 
trochemical  reaction  and  from  membrane  due  to  ohmic  loss. 
Therefore,  change  of  current  load  will  lead  to  change  of  cell 
temperature.  Fig.  9  shows  response  in  average  current  density  to 
step  change  in  cell  voltage  with  considering  about  heat  transfer. 
Comparison  with  isothermal  model,  average  current  density  of 


steady  state  at  time  t  =  1.0  s  is  increased  from  0.997  A  cm-2  to 
1 .006  A  cm-2 .  It  is  clear  from  these  result  profiles  that  increasing 
load  current  will  lead  to  increase  cell  temperature.  As  the  results, 
cell  performance  will  be  improved  due  to  the  facts  that  mem¬ 
brane  ionic  conductivity  will  be  increased  with  increasing  of  cell 
temperature  [20],  Fig.  10  shows  response  in  average  temperature 
of  membrane  to  step  change  in  cell  voltage.  The  temperature  is 
moved  smoothly  even  if  charge  double-layer  effect  is  not  con¬ 
sidered.  It  shows  that  charge  double-layer  almost  has  no  effects 
on  transient  response  of  cell  temperature.  Fig.  1 1  presents  tem¬ 
perature  distribution  on  midway  section  of  cell  for  three  time 
steps  while  cell  voltage  is  changed  from  0.6  to  0.4  V.  It  shows 
that  the  maximum  temperature  is  located  in  cathode  CL  and  the 
temperature  deviation  is  increased  with  increasing  of  current 
load. 

4.  Conclusion 

In  this  paper,  a  three-dimensional  transient  model  of  high 
temperature  PEM  fuel  cell  was  presented  to  study  the  intricate 
transient  response  to  step  change  in  cell  voltage.  Effects  of  elec¬ 
trochemical  charge  double-layer,  gas  transport,  and  temperature 
influence  on  membrane  ionic  conductivity  were  considered.  A 
single-channel  PEM  fuel  cell  with  PBI  membrane  was  adopted 
for  parametric  study.  The  contour  of  local  current  density,  oxy¬ 
gen  mole  concentration  and  electrolyte  phase  potential  were 
presented  and  discussed. 

The  results  show  that  for  the  particular  operating  conditions 
and  property  parameters  used  in  this  study,  when  charge  double¬ 
layer  effect  is  not  considered,  current  overshoot  or  undershoot 
can  be  found  with  step  change  of  cell  voltage.  Peaks  of  over- 
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shoot  and  undershoot  depend  on  air  stoichiometry  of  cathode 
side,  which  can  be  reduced  by  increasing  the  value  of  air  stoi¬ 
chiometry.  Current  density  is  moved  smoothly  and  there  is  no 
overshoot  or  undershoot,  while  charge  double-layer  effect  is  con¬ 
sidered  with  certain  value  of  specific  surface  and  differential 
capacitance  of  porous  CL.  Transient  response  of  temperature  is 
moved  smoothly  and  charge  double-layer  has  no  effects  on  it. 
The  maximum  temperature  is  located  in  cathode  CL  and  tem¬ 
perature  deviation  is  increased  with  increasing  of  current  load. 
These  transient  responses  of  PEM  fuel  cell  captured  herein  for 
the  first  time,  including  charge  double-layer  and  thermal  effects 
on  output  current,  are  expected  to  be  useful  for  design  of  power 
electronics  and  control  algorithm  for  fuel  cell  systems. 
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